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ABSTRACT 

The nature of ultracompact H ll regions (UCHRs) remains poorly determined. 
In particular, they are about an order of magnitude more common than would be 
expected if they formed around young massive stars and lasted for one dynamical 
time, around 10^ yr. We here perform three-dimensional numerical simulations 
of the expansion of an H ll region into self-gravitating, radiatively cooled gas, 
both with and without supersonic turbulent flows. In the non-turbulent case, 
we find that H ll region expansion in a collapsing core produces nearly spherical 
shells, even if the ionizing source is off-center in the core. This agrees with 
analytic models of blast waves in power-law media. In the turbulent case, we find 
that the H ll region does not disrupt the central collapsing region, but rather 
sweeps up a shell of gas in which further collapse occurs. Although this does 
not constitute triggering, as the swept-up gas would eventually have collapsed 
anyway, it does expose the collapsing regions to ionizing radiation. These objects 
can have radio flux densities consistent with unresolved UCHRs. We suggest that 
these objects, which will not all themselves form massive stars, may form the bulk 
of observed UCHRs. As the larger shell will take over 10^ years to complete its 
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evolution, this could solve the timescale problem. Our suggestion is supported 
by the ubiquitous observation of more diffuse emission from compact H ll regions 
surrounding UCHRs. 

Subject headings: H II regions 



Introduction 



When a massive star begins to emit ionizing radiation it quickly ionizes out to its initial 
Stromgren radius in the ambient gas. In typical molec ular clouds, this gas is clumpy at a ll 



scales and appe ars to be supersonically tu rbulent (e.g. iFalgarone. Puget. &: Peraultlll992l ). 
In earlier work (ILi. Mac Low, fc Abelll2004l ). we computed the initial ionization of turbulent 
gas prior to its dynamical expansion, demonstrating that the amount of mass initially ionized 
depends on the strength of the density fluctuations caused by the turbulent flow. We here 
extend that work by computing the subsequent dynamical expansion of the H ll region into 
turbulent, self-gravitating gas, driven by the overpressure of the ionized gas. 



The case of expansion into a uniform gas was first described by iKahnI (119541), and simple 



one- and two-dimensiona 



1979: Bodenheimer et al 



configurations have been considered in some de t ail (ITenorio- 



1979 



Sandford et al 



1982 



Tenorio-Tagle 



199nl : ICarcfa-Segura fc Francol 



1982 



Yorke et al. 



Tagle 



1982 



Sandford et al 



1984 



Franco et al 



1996 



Franco et al 



2007 



Arthur fc Hoarel"2006). Expansion into a u niform, magnetize d medium has been simu 



lated by lKrumholz. Stone, fc Gardinerl (120061 ). Only recently did lMellema et al.l (120061 ) com- 
pute dynamical e xpansion into a turbulent medium, though still not including self-gravity. 
Dale et al.l (120051 ) in pioneering work published the first model of dynamical expansion into 
a self-gravitating, turbulent flow, though with a model aimed more at understanding the 
global evolution of molecular clouds, as opposed to the smaller scales that we model. The 
major astronomical issue that we address here is the nature and lifetime of ultracompact 
H II regions (UCHRs). 



1.1. Ultracompact H ll Regions 



Ultracompact H ll regions (UCHR) have radii i? < 0.1 pc, and emission measures at 
ce ntimeter wavelengths of E M > 10^ cm~^ pc. Their properties have recently been reviewed 
by IChurchwelll ( 119991 . |2002| ). Their observed emission measures and sizes require that they 
be ionized by stars of type earlier than B3. If they expand at the sound speed of ionized 
gas, Cj ~ 10 km s~^, they should have lifetimes of roughly 10^ yr. Less than 1% of an OB 
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star's lifetime of a few megayears should therefore be spent within an UCHR, so the same 
fraction of OB stars should now lie within UCHRs. 



However, IWood &: Churchwelll (Il989l ) surveyed UCHRs and found numbers in our Galaxy 
consistent with over 10% of OB sta rs being surrounded by them, or equivalently, lifetimes 



> 10^ 



yr. 



Comeron fc Torral (119961 ) argued that neglecting the higher densities of massive 
stars in the molecular ring gave artificially high UCHR lifetimes. Nevertheless, they also 
derive a lifetime of 5.4 x 10^ yr, and note that even this low lifetime depends on the use 
of a high local density of ma ssive s tars d e rived from the initial mass function (IMF) of 
Humphrevs fc McElrovl Jl984jl The IScalol (119861) IMF pr edicts lower densities and longer 



lifetimes, more consistent with lWood &: Churchwelll (119891 ). 



A number of explanations have been proposed for this lifetime problem, including ther- 
mal pressure confinement in cloud cores, ram pressure confinement by in fall or bow shocks , 
champagne flows, disk evaporation, and mass- loaded stellar winds (see IChurchwelll Il999l ). 
Several of these explanations have basic problems that suggest they likely cannot explain 
the lifetime problem. 

Confinement by therm al pressure of the surrounding molecular gas requi res pressures 
of P/k ~ 10^-10^ cm-3 K JPe Free et al.lll995l : ICarcfa-Segura fc Francolll996h . At typical 
mol ecular cloud temperatures of 10-100 K, this implies densities n > 10^ cm~^. However, 
the iJeand (119021 ) mass 



M, 



-1/2 



2Gfl) 
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n 



106 



cm 



-1/2 



T 



100 K 



3/2 



(1) 



where we have assumed a mean mass per particle fi = 3.87 x 10^^^ cm~^ appropriate for 
fully molecular gas with one helium atom for every ten hydrogen nuclei. Therefore cores 
massive enough to form OB stars contain multiple Jeans masses and are thus very likely to 
be freely collapsing (e.g. iMac Low &: Klessenll2004l ). The free-fall time 
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(2) 



Typical lifetimes of > 10^ yr would thus require massive cores to last > 3ts at the hy- 
pothesized densities, rather than dynamically collapsing. Although these high pressures are 
indeed observed, they are unlikely to occur in objects with lifetimes long enough to solve the 
problem. 



A further objection was raised by lXie et al.l (119961 ). who argued t hat such high thermal 
pressures would lead to emission measures higher than those observed. (jPe Free. Goss. fc Gaume 
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19981 argue that most observed reg ions actually do agree with the predicted emission mea- 



sures, however.) IXie et al.l (119961 ) instead proposed a variation on this theme: confine- 



ment by turbulent rather than thermal pressure. However, turbulent motions decay quickly, 
with a characteristic timescale of less than a free-fall time under molecular cloud conditions 



( IStone. Ostriker. fc Gammidll998l : iMac Lowl 119991 ). Turbulent pressure would thus have to 
be continuously replenished to maintain confinement for multiple free-fall times, which would 
be difficult at such high densities and small scales. 

Another option is ram pressure confinement of UCHRs by infall of surrounding gas. 
However, this is unstable for two different physical reasons. First, the density and photon 
flux will follow different power laws in a gravitationally infalli ng reg i on ionized from within. 



so th ey can never balance each other in stable equilibrium (lYorkd Il986l : iHoUenbach et al. 



19941 ). Either the ionized region will expand, or the infall will smother the ionizing source. 
Second, the situation is Rayleigh- Taylor unstable, as this option requires the rarefied, ionized 
gas of the UCHR to support the infalling, dense gas. 



Bow shock models (jVan Buren et al.lll990l : iMac Low et al.lll99ll : [Arthur fc Hoardl2006l ) 
require hig h values of ram pressure Pram oc nv^. At n = 10^ cm~^, a velocity of ~ 10 km s~^ 
is required (jVan Buren et al.lll990l ). A star moving at such a high velocity in a straight line 
would travel a parsec over the supposed UCHR hfetime of 10^ yr, requiring a uniform- density 
region of mass > 5 x 10^ M0 for confinement. As collapse would occur on the same timescale, 
more mass would actually be required to solve the lifetime problem. Competitive accretion 
models in which massive stars orbit in the centers of newly formed groups of stars accreting 
the d ensest gas (jZinneckerl Il982l : iLarsonl Il982l : iBonnell et al.l 119971 : iBonnell. Vine, fc Bate 
2OO4J ) may address this criticism by reducing the size of the high-density region required. 



Expansion of an H ll reg ion in a density grad ient can drive supersonic champagne fiows 
down steep enough gradients ( Tenorio-Tagle 1979h. Two-dimensional rn odels first studied ex- 
pansion across sharp density discontinuities (iBodenheimer et al.lll979f). but t hen examined 
other configurations such as freely- coll apsing cloud cores JVorke et al.1ll982h. clouds with 
powe r-law density gradients in spherical (jPranco et al.lll990l ) , and cylindrical (IGarcia-Segura &: Franco 
19961 ) configurations, and exponential density gradients (lArthur &: Hoard l2006l ). Stellar 
wind combined with a champagne fiow down a continuous gradient was also modeled by 
Arthur fc Hoard (120061 ). This model does fit well the observed velocity structure in some 
cometary UCHRs such as G29.96— 0.02. However, these models face the same timescale 
problem as thermal pressure confinement models: regions dense enough to explain the ob- 
servations are gravitationally unstable, and collapse on short timescales. 

A final class of r nodels relies on mass-loaded stellar winds to reproduce the observed 



properties of UCHRs (jPyson. Williams. &: Redmarull995l : iRedman. Williams. &: Dysonlll996 
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Williams. Dyson. &: Redmanlll996l : lLizano et al.lll996l ). In these models, an expanding stellar 
wind entrains a distribution of small, self-gravitating, pressure-confined clumps that take 
substantial time to evaporate. These mo dels can reproduce many of the basic features of the 
observations including some line profiles ( Dvson et al. 1995h. core-halo, shel l ( Redman et al. 



19961 ). and cometary and bipolar shapes (IRedman. Williams, fc Dysonlll998l ). but at the cost 
of requiring an arbitrary distribution of pre-existing clumps that cannot be self-consistently 
predicted. 

In this paper we examine the expansion of H ll regions into smoothly collapsing cores 
and into turbulent, self-gravitating gas typical of a massive star-forming region. Off-center 
ionizing sources in non-turbulent cores form surprisingly round H II regions. At densities 
high enough for massive stars to form, the expanding shell driven by newly ionized gas 
quickly becomes gravitationally unstable ( jVoitlll988l : iMac Low fc Normanlll993l ). collapsing 
even more promptly than the surrounding gas. We demonstrate with numerical simulations 
that these regions of secondary collapse in the shell may be externally ionized to form objects 
with the properties of UCHRs. As the shell expands to larger sizes, new regions can form, 
extending the lifetime during which UCHRs remain visible well beyond the expansion time 
of the original H ll region. 



1.2. Summary 

In §[2] we describe our numerical methods and problem setup. In §[3] we present numerical 
results. In § H] we present analytic discussions of gravitational instability of the swept-up 
shell and the predicted radio flux from collapsing regions of the shell. We compare our results 
to observations in § O and we summarize caveats and conclusions in § O 



2. Numerical Methods 



2.1. Gas Dynamics 



We compute fully three-dimensional models of the expansion of ionization fronts in 
turbulent, gravitational l y collapsing ga s using a modified version of the code Zeus-MP vl 
( Stone fc NormanI I1992I : iNormaru |2000| ). This code uses second-order montonic advection 
(jVan Leerlll977l ). with shocks resolved using a von Neumann artificial viscosity. These al- 
gorithms are implemented in a domain-decomposed, parallel, code using the Message Pass- 
ing Inter face. The Poisson equation fo r self-gravity is computed using a F ourier transform 
method (IBurkert &: Bodenheimerlll993l ) parallelized with the FFTW library (IFrigo fc Johnson 
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20051 ). In the models of massive star- forming regions shown here, magnetic fields are ne- 



glected, as the collapsing regions are typically highly supercritical. 



2.2. Ionization 



The ionization algorithm is a parallelized version of that found in lAbel. Norman, fc Madau 



(119991 ). in which the static radiative transfer equation, 

n ■ Vly = r]u- Xvh (3) 
is solved for a monochromatic specific intensity assumed to be from a point source (/^**), 

h ■ v/f^ = (4) 

where n is a unit vector along the ray, rjy is the emission coefficient, and Xv is the absorption 
coefficient. This equation can be solved by casting rays of photons outward from the point 
source and integrating along them. We thus implement the rays in spherical coordinates 
superimposed on our existing uniform, Cartesian grid such that each ray is given two angles 
(^, (/)). We choose the number of rays such that every zone on the outermost edge of the 
computational domain receives at least one ray. We do not rotate our spherical coordinate 
system during the ru ns. This results in sma ll density artifacts appearing directly along grid 



lines, as also seen by lKrumholz et al.l (120061 ). 



Our parallelization design is rather primitive. Zeus-MP is fully domain-decomposed in 
each of the (x, z) directions. We refer to the subdomain parceled to each processor as a 
tile, and each grid point as a zone. We require that the point source be on the grid, so it will 
lie on one tile. The processor computing this source tile then sets up the spherical coordinate 
system necessary to cast rays. Each ray is then followed zone-by-zone, with each ionization 
event removing a photon from the ray, until either a tile boundary is reached or the ray is 
completely extincted. All rays that survive to tile boundaries are passed to the appropriate 
neighboring tile, where the ray walking continues. The rays at each boundary are passed as 
a batch to avoid overwhelming communications costs. 

However, this design effectively has a large serial component, because all rays must 
be propagated across the source tile before they can be computed on other processors. As 
a result, parallelization is limited to rou ghly 32 processors. M ore efficient parallelization 



methods, such as adaptive ray tracing (lAbel fc Wandeltl |2002| ). are clearly necessary for 
larger scale calculations. 



Ionization is followed with a passive tracer field X that is advected as an intensive 
quantity such as temperature, rather than an extensive quantity such as density. In practice. 
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this is implemented by advecting the quantity pX and then dividing by the updated density p 
after the advection step to recover the tracer field. The ionization is initialized to T = 0; when 
zones are computed to be ionized during the ray tracing, they are set to T = 1. Intermediate 
ionization values occur due to numerical diffusion during advection. However, this procedure 
assumes that recombination happens on timescales shorter than the dynamical timescale, so 
no partially recombined regions will be traced. 



2.3. Heating and Cooling 



We explicitly incorporate a very simple heating and cooling model, allowing ionization X 
to control the heating rate. We solve an implicit equation for the internal energy at timestep 
m + 1 

= - l)e-+iV • V - Cin^Ale^+i) + nT{X)l (5) 

where At is the timestep, 7 the adiabatic index, v the velocity, n the number density, A 
the temperature dependent cooling rate, V the ionization dependent heating rate, and C a 
constant that we include for numerical convenience as we discuss below. To solve the implicit 
equation, we use a Newton- Raphson method with a binary search when that fails to converge 
(IPress et al.lll992h . 

Our models cannot resolve the cooling lengthscales of 10^^ — 10^'^ cm typical of dense gas, 
nor will it be computationally feasible to do so in the near future for models encompassing 
an entire compact H ll region. Therefore we only implement a simple approximation to 
the cooling function A(T), b ased on a modified version of the radiative losses table from 



Flash v2.3 (jPrvxell et al.ll2000[) . The values of that table at T < 10'* K are taken from 
Dalgarno &: McCrav (|l972), ass umin g an ionizatiori fraction of 10~^, and at T > 10^ K from 



Raymond. Cox, fc SmithI (119761 ) and ISarazinI (Il986l ) . This captures the basic behavior that 
neutral gas quickly cools to a uniform temperature, although the ionization fraction assumed 
is more appropriate for diffuse gas than molecular cloud gas, the radiative cooling rate for 
ionized gas depends on the ion and electron densities and rii rather than n, and we are using 
a temperature-dependent equilibrium cooling curve rather than following the non-equilibrium 
behavior of the ionized (IBenjamin. Benson, fc Coxll200ll ) and neutral (IGlover fc Mac Low 
20071 ) gas. Our main purpose here is to capture the qualitative dynamical behavior of the 
expanding ionized region, however, and for that this approximation is sufficient. 



Following a similar philosophy, our heating function is empirically chosen to fix the 
equilibrium temperature in the neutral and ionized gas at the initial density hq by setting it 
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to 

r = A(T,)J + A(TO(l-T), (6) 

where Tj = 10^ K is chosen to approximate the ionized gas temperature, and T„ = 100 K 
is chosen for the neutral gas temperature in the photodissociation region surrounding the 
ionized region. (Note that we are taking the magnitude of A, but not its dimensionahty, as 
we are still assuming that the heating term is nT as appropriate for photoionization heating.) 
The moderately high neutral temperature also serves to keep th e Jeans mass higher, making 
it easier to resolve gravitational collapse (iTruelove et al.l 119971 ). To resolve the heating and 
cooling, we must add a condition to our timestep requiring 



At < Cce{t)/{de/dt) 
where we set Cc = 0.3, and {de/dt)rad = C(^^^ ~ ^^)- 



(7) 



However, for molecular cloud densities and the temperatures we have chosen, this 
timestep can be as much as six orders of magnitude shorter than the dynamical Courant 
timestep At < CAx/ max(v,Cs), where C ~ 0.5 is the Courant number. As our primary 
purpose, again, is to qualitatively capture the behavior of the expanding high-pressure H ll 
region, we choose to set ( = 10~^ in equation ([5]). This has the effect of slowing the ap- 
proach to the equilibrium te mperature, and a l so of increasing the thickness of the cooling 
region behind strong shocks. lArthur fc Hoard (120061 ) used a similar approximation in their 
two-dimensional models. 

We have run tests of the expansion of spherical H ll regions (model A) to examine this 
approximation. The radius of the H ll region is shown in Figure [1], which shows that, for 
the densities chosen in model A, small values of ( begin to significantly impact the behavior 
of the H II region. (We terminated the C = 0.1 run early because of the long computation 
time — several weeks — required even for this modest resolution.) The chosen value or larger 
results in less than 10% errors in the radius over time, which we took to be sufficiently 
accurate for our purposes. This is s hown in Figu re [H which compares the numerical solution 



to the usual analytic solution (e.g. ISpitzerlll978l ) for values of C = 10 ^, 0.01, and 0.1. (The 



actual radius measured in the numerical solution is one zone less than the first zone with 
color greater than 0.5 along the 

Most of the error is at early times when the lower temperatures due to the slowed 
heating rate at low values of ( slows the expansion of the H ll region. Figure [2] shows 
the central temperature as a function of time for these three models, showing the slow 



approach to equ i. 



ibrium of the low C model. This appears to contradict the speculation of 



Krumholz et al.l (120061 ) that the error is primarily due to enhanced cooling at the ionization 



front. We note that model A uses lower density and ionizing luminosity than our production 
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models. As the heating and coohng timescales depend on the density, this means longer 
thermal timescales in model A than in the other models. The actual values we chose for 
our production models yield a far faster approach to equilibrium temperature, though, as 
shown by the dotted curve in Figure [21 which follows a model using the same density and 
ionizing luminosity as models E-J. This model is also shown in Figure [H Thus, our choice 
of ( = 10~^ appears appropriate for our production models, though marginal for our test 
case model A. 



2.4. Problem Setup 



Our goal is to simulate the evolution of H ll regions produced by single stars in realistic 
density distributions produced by gravitational collapse, taking into account off-center stellar 
positions. We describe two different types of models, with parameters given in Table 1. Both 
are computed in cubes with periodic boundary conditions on all sides. 

We first examine off-center stars in otherwise smoothly collapsing cores, with no further 
structure. We choose to study off-center stars in the smooth case because the centered 
case has been previously treated (lYorke et al.lll982l : iGarcia-Segura fc Francdll996l ). with the 
natural result that centered stars yield spherical H ll regions unlike most observed UCHRs. 
A similar argument leads us to allow the stars to begin ionization off-center from the density 
peak in the turbulent case. We begin with gas uniformly distributed on the grid, with 
a 1% density perturbation in a central sphere to ensure the core collapses in the center. 
The dynamically collapsing cores in these models develop density profiles with flattened 
centers, as shown in Figure [31 and the stars are placed at varying distances from their centers. 
The peak density is centered in the box. Note that the inner flattening occurs at radii of 
only a few zones, and is therefore likely numerical rather than physical. 

We then examine turbulent model s. They ar e set u p by first uniformly driving turbu- 
lence using the algorithm described by [Mac Lowl (Il999l ) , with driving parameters given in 
the caption to Table 1. This gives an rms velocity of Vrms = 1-5 km s~^, and an rms Mach 
number of Airms = 2.4. 

We drive for 0.7 Myr, several times the crossing time of the largest eddies with size L/2, 
sufficient to establish a steady-state flow, before turning on self-gravity and allowing collapse 
to begin. The Jeans number Nj = M/Mj of our models is given in Table 1, where M is the 
mass contained in the computational domain, and the Jeans mass is given by equation ([T|), 
so that 

0.2kms-V \PoJ V2pc 



Nj = 267 Mr, 
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The density scaling po = 1-928 x 10~^^ g cm~^. We choose large values of Nj appropriate for 
regions undergoing massive star formation. Our computational domain contains sufficient 
mass that it also exceeds the turbulent Jeans mass Mj^t = {vrms/ CsY^j by at least an order 
of magnitude, so collapse occurs quickly when gravity is turned on. 

We turn on ionizing rad i ation once collapse has proceeded to the limit of our resolution 



following the iTruelove et al.l (119971 ) criterion that the local Jeans length be resolved by at 
least four zones. For typical models, the actual collapsing core at that point only contains 
several Jeans masses, so our models are of relatively small stars of 3-4 Mq. We intend them 
primarily as experiments to reveal the basic behavior of an ionization front under these 
circumstances, rather than as models of specific objects. 

The largest of these runs, model H, ran for 54 days on 32 processors. Because of the lim- 
itations of our parallelization algorithm, we cannot gain substantial additional performance 
by going to higher processor number. 



3. Results 

3.1. Smooth Collapse 

We first consider a model of an ionizing region off-center in a collapse in a uniform flow. 
Here we can clearly follow the morphology o f the champagne flow for med in a gravitationally 



collapsing region. Unlike previous work (e.g. iFranco et al.lll990l . 120071 ). we follow dynamically 



both the gravitational collapse and the expansion of the ionized gas in three dimensions. 
We turn on the ionizing source at various positions off the center of the collapsing core 
after collapse has proceeded to form a spherical core surrounded by an envelope with radial 
density dependence p oc (Fig. [3^). 

In Figure H] we show the time development of the H II region resulting from turning 
on an ionizing source 0.22 pc diagonally from the center of the core (model D). As can be 
seen from Figure [3^, this lies within the power-law envelope. Varying the position of the 
ionizing source from the inner edge of the envelope at 0.125 pc to the middle of the envelope 
at 0.22 pc (models B, C, and D) makes little difference to the conclusions drawn here, as 
shown in Figure [5l 

The expansion of the H ll region drives a strong shock into the surrounding gas, sweeping 
up a thin, dense shell of neutral gas that traps the H ll region for the 10^ yr duration of 
these runs. In the smooth core model shown here, the neutral sound speed is 0.2 k r n s"\ 



corresponding to gas with a temperature of order 10 K. Recent work by lFranco et al.l (120071 ) 
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uses an effective sound speed more than an order of magnitude liiglier as an approximation 
to turbulent motions in order to maintain liydrostatic equilibrium in their cores, rather than 
following the collapse as we do. As a result they do not see dense shell formation. 



Contrary to much early speculation (e.g.. lFey et al.lll992l : iGaume et al.lll995l ). the strat- 
ification of the spherical envelope produces a nearly spherical H ll region shape during the 
period that we simulate, rather than a cometary shape. Th at this migh t occu r was actu- 
ally fii;st_jhown_in a dif ferent context by the analytic work of iKorycanskyI (119921 ). who used 
the iKompaneetd (Il960l ) approximation to compute the shape of off-center blast waves in 
spherical, power-law stratified, density distributions. Confinement in conical regions only 
happens for density power laws steeper than -4. Blast waves expanding into distributions 
with power laws in the range -4/3 to -8/ 3 ultimately open out and wrap completely around 
the central core. [Arthur fc Hoard (120061 ) find a similar spherical shape to our own result in 
their model F, which al s o incl udes a stellar wind. This result has been more comprehensively 
investigated by lArthun (120071 ) who shows that nearly spherical H ll regions are expected in 
a wide varie ty of power-law dens ity distributions, so long as the H ll region does not become 
unbounded. iFranco et al.l (120071 ) found cometary shapes, but it was most likely because they 
computed two-dimensional models assuming slab symmetry: effectively a line ionizing source 
near a cylinder. The lack of the possibilit y of expanding in the third dimensiori tends to give 
narrower blowouts (this was also seen by lMac Low. McCray. fc Normanlll989l in a different 
situation) . 

In the models presented here, the innermost region of the collapsing core is too dense 
to be ionized. It continues to collapse despite the ionization of its envelope, as shown 
in Figure [61 Indeed, the impact of the photoionization-driven shell appears if anything to 
accelerate collapse, as our model B, with the ionizing source closest to the center of the cloud, 
has a higher peak density over time. However, note that model D, with source somewhat 
farther from the peak than model C, still ends up with slightly higher density. This may 
just be the effect of the ray-tracing artifacts which impact the core in model C, and not in 
model D. A nearby ionizing source will in general have great difficulty preventing collapse in 
a region that has already reached moderately high density because the recombination rate 
oc 'n? while the ionization rate is only cx n, so as collapse proceeds, more and more photons 
are required to maintain the ionization of the same region. We discuss this quantitatively 
below, in § lOl 

Although the shape of the shock front is close to spherical, the intensity of emission 
from the neutral shell and the ionized interior varies around the shell, as the shell is higher 
density where it hes nearest to the center of the core (Fig. H]). We do not explicitly show 
simulated observations of the ionized gas, as artifacts in the boundary layer ionized density, 
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caused by recent ionization of high-density neutral zones, dominate the emission. Instead, 
we show in Figure [7] the column density J ndi, which is dominated by the neutral density. 

The highest column density point is at the dense center of the gravitationally collapsing 
core. Even before this forms a star, as it reaches higher densities and pressures, the gas 
ionized on its surface also reaches higher densities. This ionized gas may have sufficiently 
high emission measure to appear as a U CHR, surrounded b y a more diffuse compact H ll 
region, the usual configuration observed (IKim fc Kod 120011 ). However, it will have a brief 
lifetime, of order the free-fall time of 10^ yr or less (see eq. [2]). In the next subsection we 
consider the development of a more complex fiow. 



3.2. Turbulent Flow 

To understand the development of H ll regions in a more realistic environment, we 
examine their formation in cores that are self-consistently collapsing from the supersonic 
turbulent fiow characteristic of a molecular cloud. We turn on ionization near the highest 
density point in the simulation, at the time that the density reaches the Jeans resolution 
limit. Figure [8] shows that only our highest resolution model remains formally well resolved 
after collapse begins, but that the behavior of the peak density is reasonably well resolved. 
To follow the expansion of the H ll region for the longest time possible, we take advantage 
of the periodic boundary conditions to shift the highest density point to the center of the 
cube before beginning ionization. Once again, the density structure of the cores takes on a 
radial density dependence p oc (Fig. [3b). 

In Figures M and [TU] we show the expansion of the resulting H II region in our highest 
resolution model H, with the ionizing source only 0.05 pc from the center of the collapsing 
core, in the inner portion of the power- law envelope. In Figure [TO] we compare the morphology 
of models E (128^ zones) and H (256^ zones). The primary differences are in the thickness of 
the shell, which is not fully resolved even at the higher resolution, and in the wavelengths of 
instability resolved in the shell, again probably not fully resolved. However, the qualitative 
results that we discuss are ones that the two models agree on. At the lower resolution 
we varied the position of the ionizing source in models E, F, and G, but, as shown in 
Figure [Til this again made little qualitative difference. These different models run with the 
same background density field do give some idea of how sensitive our results are to small 
perturbations. 

The expanding blast wave driven by the ionized gas encounters strong ambient density 
fiuctuations produced by the background supersonic turbulent fiow. These fiuctuations both 
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directly shape the shell and seed instabilities in the expanding shell. Four instabilities are 
evident. In regions where the shell is expanding into a sufficiently low density region to be- 
gin accelerating, Rayleigh- Taylor instabilities occur. (Strictly speaking, as these are driven 
by acceleration rather than gravity, they should be denoted Richtmeyer-Meshkov instabili- 
ties, as is commonly done in the fluid dynamics comr nunity.) Second , a thin, decelerating, 
pressure-driven, shock-bounded shell is subject to the Vishniac ( 1983h instability. Our res - 
olution is probably insufficient to capture its saturated state (IMac Low fc NormanI Il993l ). 
but some clumping se en in the shell will be caused by even an underresolved manifestation 
of the instability (cf. iMac Low et al.lll989f). Thi r d, the introduction of ionizing radiation 



impingi ng on a decelerating shel l drive s iGiulianil (119791 ) instabilities, as numerically mod- 
eled by iGarcia-Segura fc Francol (Il996l ). Finally, and perhaps most interestingly, the shell 
itself quickly accumulates enough mass at high enough densities t o become gravitationally 
unstable, something pred icted for many years by analytic models ( lElme green &: Ladal 119771 : 
Vishniacl[l983l : IVoitlll988h . 



Although the ionized gas readily expands into low density regions, it does not appear to 
disrupt the collapsing core. Figure M shows no perturbation to the increasing density caused 
by collapse when ionization turns on at t = 0.83 Myr. This conclusion is supported by the 
analytic computation presented below (§ 14. 2p . The numerical result must be approached 
with some caution, though, as the Jeans length in the center of the core is resolved by less 
than four zones during the period after ionization begins. 

As gravitational instability in the shell sets in, collapse be gins at multiple points. To 
charac terize this behavior, we use d the clumpfind algorithm of IWilliams. de Geus. fc Blitz 
(119941 ). optimized as described by lKlessen. Heitsch. fc Mac Low to flnd clumps with 



density at least 10% of the peak density of 5.7 x 10~^^ g cm~^ in the cube at the last timestep 
of model H. Contours of 5% of the peak density were used. This captures only regions that 
are undergoing gravitational collapse, as can be seen from Figure [HI which shows that prior 
to turning on ionization and gravity, peak densities are approximately 1% of the flnal peak 
value. Even accounting for shock compression of those peak densities, we show below in 
§ 14.11 that the peak shell density absent gravitational instability will be only 7% of the peak 
density. 

Aside from the original collapsing core, we flnd two other collapsing regions with masses 
of a few solar masses, of the same order of magnitude as the original core. Again, we 
must emphasize that this is a qualitative result, as these regions have collapsed beyond 
the Jeans limit. However, they are sufficiently well separated by resolved gas to rule out 
spurious fragmentation in the shell as causing them. Whether we have fully resolved their 
internal structure is another matter: we actually found 10 clumps distributed among the 
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three regions, but do not consider this sub-clumping to be well-resolved. We discuss the 
instability of the shell further in the next section. 

In Figure [T2] we show the column density distribution of the gas in our models. As 
in Figure [TJ this is dominated by neutral molecular gas, while the ionized gas shows up 
as low column density cavities. There are several points to note about the column density 
distribution. The material swept up in the neutral shell surrounding the ionized region is dis- 
tinguishable from the background turbulent gas morphologically. The shell is thinner, with 
higher contrast and more small-scale structure. This occurs primarily because the back- 
ground turbulence is at a substantially lower Mach number than the shell, so the structures 
formed in it are thicker and have lower density contrast. The regions of secondary collapse 
can also be picked out in the column density projections. In particular, in the yz projection, 
the regions to the left are clearly separated from the primary core in the center of the image, 
although they are projected almost on top of each other in the xz projection. 



4. Analytic Considerations 
4.1. Gravitational Instability 



We can understand the observed gravitational collapse by examining the c ollaps e cri- 
terion for expanding sh e lls. T his was first computed by lElmegreen fc Ladal (119771 ) and 
Elme green &: Elmegreenl (119781 ). who assumed that turbulent velocit ie s in t he sh ell would 
be of the same order as the expansion velocity. lOstriker fc Cowid (jl98ll ) and IVishniac 
( 1l983l ) assumed, on the other hand, that the turbulent velocities in the shell would only 



be transonic. Thi s latter assumption was sup ported by analytic work by IVoitI (119881 ) and 
nu merical work by Mac Low fc NormanI ( 19931 ). We here use the formulation of this criterion 



by lMcCray fc Kafatod (Il987l ). who demonstrated that collapse will occur in a spherical shell 
of radius i?^, expansion velocity ambient density p, and sound speed Cg if the gravitational 
potential energy in a segment of shell exceeds its pressure and kinetic energy of expansion, 
a condition that can be expressed as 



T = {).QlGpRl/{VsCs) > 1. 



(9) 



In our case the shell expands into a medium of varying ambient density p. However, this 
criterion is derived locally for each patch of shell, so approximating the ambient density with 
the local value is a reasonable assumption. (Although the ambient density varies radially as 
well as angularly, most of the mass in any patch of shell accumulates during passage through 
high density regions near the current radius.) 
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Our models show shell collapse away from the primary core in the turbulent case, but 
not in the smooth case. The turbulent model H has a maximum ambient density of p ~ 
5 X 10~^^ g cm~'^ (Fig- ED, a shell radius Rg — 0.25 pc and velocity — 5 km s~^, and 
sound speed = c„ = 0.63 km s~^, giving T ~ 4. Parts of the shell hitting dense filaments 
in the background turbulent flow become gravitationally unstable and collapse, while parts 
expanding into lower ambient densities remain stable. The smooth model B, on the other 
hand, has an average ambient density p ~ 6 x 10~^^ g cm~^ with less fluctuation in density 
from the average than in the turbulent case. Its radius at the end of the run is Rg ~ 0.5 pc, 
and sound speed c„ = 0.2 km s~^, giving a value of T ~ 0.06, in agreement with the lack of 
collapse found in the model. 

The Jeans mass in unstable portions of the shocked shell of model H can be estimated by 
taking the peak density in the ambient gas, and assuming that it is hit by an isothermal shock 
of Mach number A4 = Vg/cn ~ 8, so that the shell density ps = Ad'^p ^ 4 x 10"^'' g cm~^, 
or, in terms of number density, Ug — 10^ cm~^. Substituting into equation ([1]), we find 
Mj ~ 2 Mq, agreeing within better than 50% with the masses of the observed regions of 
secondary collapse in model H. 



4.2. Trapping Radius 

At what radius does an external ionization front get trapped in a collapsing core? We 
can calculate the initial standoff distance of an ionization front impinging on a collapsing 
core. If Tg is small, the core is quickly ionized, while if approaches the distance to the 
ionizing source d, the ionization front will have little influence on the core. 

The freely collapsing objects formed by gravitational instability in our models have 
density profiles, while a core with constant central accretion rate M = —ATrR'^pnQVg, w here 



the f ree-fall velocity = — (2GM/i?)^/^, has a shallower r '^/^ density profile (e.g. I Shu 



1992|). 



We work in the frame of reference of the ionizing source, with distance from the source 
given by r. Along the line connecting the source with the center of the core, R = d — r. 
We take the number density n{r) = no{R/ Ro)~^ , where uq and Rq characterize the mass 
and size of the core, and 77 = 3/2 or 2. Assuming the photoionized gas to be fully ionized, 
ionizing photons are absorbed at a rate 

— = -A-nr'^nirYa (10) 
dr 



(e.g. ISpitzeii Il978l ). where a is the hydrogen recombination coefficient to the second level. 
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Substituting for density and integrating both sides, we find the stellar ionizing photon lumi- 
nosity 

= AimlR^^a I ' drr'^{d - r)-^^. (11) 
Jo 

The integral can be evaluated by substitution of variables for = 3/2 or 2. 

We can place the solutions in dimensionless form. Define the fractional standoff distance 
of the ionization front from the center of the core along the line of sight to the star ^ = 
{d — rs)/d. Small values of ^ mean highly ionized cores, while large values suggest only the 
core surface is ionized. Then we can define a dimensionless ionizing photon luminosity 



Note that for rj = 3/2, this dimensionless luminosity is independent of the actual standoff 
distance of the core d. Then, the solution for cores with density power-law r/ = 3/2 is 

A3/2) = ^-lne-^ + ^, (13) 

and for 77 = 2, it is 

A2)^4 + i-i+3^. (14) 

If we substitute typical values, including ionization coefficient a = 2 x 10^^^ cm^ s^^, we 
find 



\ f no f R{ 



and 



^3/2) ^ 1.4 ( ^ J j ^ (15) 

£(2) = (1.4 X 10^) ( -4^) ( ^) ( "° \ " ( " . (16) 

The value of jC{ri) determines whether a collapsing core is promptly ionized or remains mostly 
neutral. The higher value of jC{2) for the same parameters shows that the corresponding 
cores are more easily ionized, as expected. In Figure [T3] we show the relation between the 
fractional ionization standoff distance ^ and the dimensionless luminosity jC{ri) for 77 = 3/2 
and 2. For plausible parameters, ^ is sufficiently large to suggest that these cores could form 
high-emission measure sources for a free-fall time. 

We note that observed molecular cloud core s have power-law envelopes with flat interiors 



[e.g. lAndre. Ward-Thompson, fc Barsonyll2000l ). If the calculated stand-off distance for such 
a core becomes less than the radius of the const ant- density central region, we expect the core 
to be quickly ionized. 
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4.3. Expected Radio Flux 



We can use our knowledge of the trapping radius to make an estimate of the expected 
radio flux density fr om the ionize d portion of a collapsing core, for comparison with observed 
UCHRs. We follow iGarayl (119871 ). who assume that the flow of ionized gas from a dens e core 
embedded in an H ll region follows an rj = 2 power-law density profile. iHenneyl (120031 ) show 
analytically that this is likely to be a good approximation f or many radia l scale lengths away 
from a core embedded in the wall of an H ll region, while [Franco et al.l (120001 ) demonstrate 
that a number of observed UCHRs have radio spectral indexes consistent with power-law 
profiles of ionized density. 



The radio f 



ux den sity from an unresolved source depends on whether it is optically thin 



or thick. iMoranI (119831) gives the op tically thin flux density at frequency u, using power law 
opacity (IMezger fc Henderson! 119671 ). as 



= (3.4 mJy) ( 



1 GHz 



-0.1 



V 



10^'' cm 



10^ K 



-0.35 



D 



1 kpc 



;i7) 



where is the electron temperature, D is the distance, and V = / nldV is the volume 
emission measure, computed from the electron number density He and the volume of ionized 
gas V. The volume emission measure for an object with ionized number density n(r) = 
no(r/ro)~^ for r > ri can be shown to be V = AirnQrf. In our case, the objects in question 
are density enhancements in the shell. The one-sided nature of the ionization will make 
about a factor of two difference in the volume emission measure, which we param eterize 
with a factor tp. Combining this with equation (|T7|) . we can write the flux density as (I Gar ay 
1987h 



= (5.4 mJy) ( ^ 



5 X 10^ cm' 



1015 



cm 



1 GHz 



-0.1 



T 



10^ K 



-0.35 



D 



1 kpc 



Using equation (|T6l) and Figure [13] we can determine the ionization radius = d{l —C,)- We 
can then find the density at that radius = nQ{ri/rQ)~'^ to derive the flux density expected 
for any given core parameters, luminosity, and distance. 



5. Comparisons to Observations 



Our models give insight into the behavior of smaller regions of massive star forma- 
tion, where man y observed UCHRs r e side, partic ularly those picked up in broad surveys 
such as those of IWood fc Churchwelll Jl989h and Eurtz^t^ (Il994l ). These UCHRs are 
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nearl y always observed t o have associated extended emission (iGaray et al.lll993l : iKurtz et al. 
I999I: iKim fc Kool I2OOI : Isik et aP boosl ) consistent with a larger underlying H ll region. 
Kim fc Kool (I2OOII ) suggested that this was simply due to the hierarchical structure of molec- 
ular clouds, but the lifetimes for the densest r egions in their scenario were derived from the 
photoevaporation timescale (jWhitworthl Il979l ) . which can be more than an order of magni- 
tude longer than the Jeans collapse timescale for these dense regions. Instead, the models 
shown here suggest that the o bserved UCHRs are short-lived regions of collapse in a longer- 
lived, larger, expanding shell. iBik et al.l (120051 ) found that many UCHRs do not themselves 
contain OB stars, but have nearby embedded clusters with OB members. As these clusters 
have ages of order 10® yr, they suggest that the cluster age, rather than the age of individual 
UCHRs, is the relevant timescale, in agreement with the model presented here. 



Unresolved UCHRs in the survey of IWood &: Churchwelll (119891 ) have integrated flux 
densities at 2 and 6 cm in the range of 2-1700 mJy. We can compare the flux density ex- 
pected for cores in our model to confirm whether it is consistent with the observed range. 
Unfortunately, the lack of a detailed treatment of the interface between neutral and ionized 
gas inhibits our ability to directly measure the ionized ma terial flowing off density condensa- 
tions in the simulated shell as predicted by iHenneyl (120031 ). Nevertheless, we can make some 
order of magnitude estimates of the properties of the clumps. 



In § 13.21 we showed that regions of secondary collapse typically had masses of a few solar 
masses, which we'll take for the example to be M = 2 Mq. To set scales, we choose the back- 
ground density of the H ll region at the end of model H of po — 2 x 10"^*^ g cm~^ as the outer 



edge of the r ^ region. Then tq = (M/47rpc 



a/3 



0.08 pc. If we take the mean mass per par- 



ticle for molecular gas with 10% He nuclei /i = 3.39 x 10^^^ cm~'^, then no = 5.9 x 10^ cm~^. 
A typical distance to the central star is (i ~ 0.2 pc, and in this model the star's ionizing lumi- 
nosity is = 10^^ s^^ Substituting into equation (ITB]) . we find C{2) = 4 x lO^'^. Examining 
Figure [T^ we find ^ ~ 0.86, and thus = d{l — ^) = 0.028 pc. The density at that radius 
is then rii = no(rj/ro)~^ = 4.8 x 10^ cm~^. We can then use equation (|T8|) to derive the flux 
density from this object, = (970 mJy)(D/5 kpc)-2(z//15 GHz)-°-i(T/10^ K)-o-35(^/0.5). 
For typical values, this is consistent wit h the flux densities of even bright unresolved objects 
observed by lWood fc Churchwell Jl989h . 



Kurtz et al.l (119941 ) found a direct correlation between size and density of spherical 



and unresolved UCHRs. They suggested that expanding regions would behave in this way. 
However, collapse offers an alternative explanation for th e observed correlatio n . Ob serva- 



tions of massive star formation r egions such as Sag B2 (IDe Free et al.l Il996l . Il998l ). W3 



dTieftrunk et~aDll997h . and W49 JPe Free et al.lll997l . [200oh show multiple UCHRs closely 
spaced. These regions are far larger than the ones simulated here. Further work will need 
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to be done to understand whether these larger regions generate a single expanding shell as 
here, or many smaller isolated ones. However, collapse in a single expanding shell would 
again address the lifetime question. 

In simple models based on static spherical cores, the calculated number of ionizing 
photons required for observed UCHRs have typically been found to be factors of 3-10 higher 
than can be provided by external ionizing sources. However, collapsing structures in a sheet 
could conceivably produce the same emission measures with fewer photons. This problem 
requires further investigation. 

Most high-resolution observations of H ll regions to date observe either emission from 
the ionized gas, proportional to the square of the ionized gas density, extinction of the ionized 
gas by the neutral gas, or near-infrared emission from small dust grains or PAHs, presumably 
embedded in the neutral gas or the densest ionized gas. The last is determined by a combi- 
nation of radiative heating from the central star and local dust density. Even sub-millimeter 
dust emission, which traces column density in cold clouds, will be strongly perturbed by ra- 
diative heating near the surface of the H ll region. The pure column density map shown here 
in Figure [12] is suggestive of the results expected from these different observational methods, 
but lacking any treatment of radiative heating or direct measure of the ionized density, does 
not represent a direct simulation of any of them. 

That caveat being stated, the shape of the cavity is certainly strongly reminiscent of the 
observations of larger H ll regions such as th ose by the GLIMPSE survey with the Spitzer 



Space Telescope (e.g. IChurchwell et al.ll2004l ). In particular, we reproduce the filamentary 
nature of the dense gas, and the fingerlike protrusions in from the wall of the shell, as well as, 
naturally, the appearance of an internal cavity. We note that, although we see high column 
density filaments, the shell is dense in every direction, with the filaments simply being the 
highest column density regions. 



6. Conclusions 



We have modeled the dynamical growth of young H ll regions in turbulent, self-gravitating, 
molecular clouds. Our assume d initia l conditions reflect the star formation paradigm re- 
viewed by iMac Low fc KlessenI (120041 ). This suggests that massive star formation occurs 
in cloud cores collapsing from a turbulent fl ow that are far from hydrosta tic equilibrium, 
as already proposed in the classic review by IShu. Adams, fc Lizand (119871 ). The collapse 
produces the observed high pressures, so they are transient, lasting only on the order of a 
free-fall time. The medium from which the collapse occurs has large density fluctuations 
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produced by a supersonic turbulent flow, so that the H ll region resulting from massive star 
formation expands into an inhomogeneous medium. 

Before giving the conclusions we draw from the simulations, we summarize the strengths 
and limitations of these models. They are fully three-dimensional models including driven, 
hypersonic turbulence, that use a ray-tracing algorithm to follow direct ionization from a 
single ionizing star. Although our highest resolution models have 256^ zones, and resolve the 
Jeans length with four zones everywhere except the very centers of a few collapsing regions 
at late times, they do not fully resolve the dense, swept-up shell behind the expanding shock 
front, nor the denser ionized gas layer coming off of collapsing cores. We use periodic bound- 
ary conditions, so large-scale gradients in the molecular gas are not modelled. Scattering of 
ionizing radiation is not accounted for, so that shadows are sharper than they should be. 
Magnetic fields are not included in these models, though they might modify the structure 
of the dense shell. They are unlikely to be strong enough to resist gravitational collapse in 
the massive star-forming region we are modeling, however. We have artificially extended 
the heating and cooling timescales, reducing the pressure in the ionized gas at early times, 
and thus making as much as a 10-20% error in the radii of the resulting H ll regions. Ab- 
sorption of ionizing radiation by dust is also n ot included, which can reduce the available 



ionizing luminosity by an order of magnitude (jPetrosian. Silk, fc Fieldl 1 19721 : [Arthur et al. 



2004 ). though probably in a uniform fashion that will not impact ou r qualitative conclusions. 
Finally, we have not included stellar winds ( lArthur fc Hoard l2006l ). which will drive faster 
expansion and more accumulation of mass in the swept-up shell, tending to counteract the 
effect of dust. 

Despite their limitations, our models give insight into the nature of UCHRs. We 
demonstrate that the shape of an H ll region expanding off-center in a core with an r^^ 
p ower-law den s ity di stribution is roughly spherical, confirming the general analytic results 
of iKorycanskyI (119921 ). This suggests that simple champagne flow mode ls may have difficulty 
prod uci ng the parabolic shape s characteristic of cometary UCHRs (IWood &: Churchwell 
19891 ). lArthur fc Hoard (120061 ) were able to form the parabolic shapes observed only in 
a planar density gradient with a stellar wind confining the H ll region. They also find a 
spherical shape in a spherical core, though (their model F). 

The expanding shell of the H ll region takes roughly 10^ yr to reach a radius of a 
parsec and break out of its parent molecular cloud. During that time, self-gravitating cores 
repeatedly collapse in it as it interacts with pre-existing turbulent density fluctuations. While 
each core collapses in only ~ 10^ yr, they form repeatedly over the lifetime of the shell, and 
can be externally ionized to emit radio fluxes comparable to the observations (§ 14. 3p . Since 
these cores will often be too low mass to form OB stars, this scenario might explain both the 



apparent lifetime of UCH Rs of > 10^ yr (IWood fc Churchwelll Il989l ) and their association 
with compact H ll regions ( Kim fc Koo 2001 : Bik et aL 20051 ) . 



Work remains to be done to demonstrate the viabihty of this scenario. Most importantly, 
the expected velocity structure and detailed emissivity structure for externally ionized cores 
in an expanding shell need to be computed. The relative importance of externally ionized 
cores in an expanding shell that might account for unresolved sources, and stars orbiting in 
a cluster potential that could produce cometary UCHRs also needs further consideration. 
These two scenarios make distinctly different predictions for the direction of ionized gas flow, 
and the relationship between ionized and molecular gas. 
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Table 1. Model Parameters 
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0.125 
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0.125 
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3 


1390 


0.125 


0.125 


0.125 
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0.05 
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256 
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10 


0.63 


100 
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0.05 
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64 
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10 


0.63 


100 
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0.05 









Note. — All models have ionized sound speed q = 10 km s~^, and heating and cooling 
reduction factor C, = 10~^. All turbulent models have driving luminosity E — 1.875 x 
2q33 gj,g g-i driving wavelengths —1-2. 

^Size of cubical computational domain in parsecs. 

'^Number of zones on a side of the computational domain. 

'^Initial conditions: either smooth (U) or turbulent (T), driven with driving wavenumbers 
kn =1-2 and driving energy input E = 1.875 x 10^^ erg s~^. 

'^Models with gravity are indicated by Y, others by N. 

^Scaled ionizing photon luminosity 5'*/(10^^ s~^). 

^Neutral sound speed in km s~^. 

^Average density p scaled by po — 1.928 x 10~^^ g cm~^. 
''Number of Jeans masses in computational domain. 

^Distance in parsecs of ionizing source (in negative direction) from peak density along 
X-axis, and correspondingly along y- and 2;-axes. Shifts are 8 zones in 128^ models. 
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Fig. 1. — Radius of ionized region in homogeneous medium over time from the theoretical 
prediction {solid line), and our model A at a resolution of 128'^ zones with values of the 
heating and cooling reduction factor ^ = 0.1 (squares), 0.01 {diamonds), and 10~^ {asterisks). 
Computation time linearly increases with ( so long as the cooling time step dominates. A 
model with the same density and temperature as models E-J is also shown crosses, along 
with its theoretical prediction dashed line. 
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Fig. 2. — Central temperature of ionized region in homogeneous medium over time from our 
model A with values of the heating and cooling reduction factor ( = 0.1 (solid), 0.01 (dotted), 
and 10^^ (dashed). The figure shows the slower approach to equilibrium temperature of the 
lower ( model that we use. However, as the thermal timescale depends directly on the density, 
we also show with the thin dashed line a model with the same density and temperature as 
models El-J, which approaches equilibrium in far less than a dynamical time. (Note that the 
dips at late times are caused by the dynamics within the H ll region driven by the rarefaction 
wave that opens behind the initial shock front.) 
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Fig. 3. — Density as a function of radius of the collapsing core immediately before ionization 
is turned on in (a) a uniform medium along the x-axis (128^ model B), and (&) a turbulent 
medium (256^ model H) along the x-axis {solid), y-axis (dashed), and z-axis (dash-dotted). 
The thin dotted line shows the power law in each case. 
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Fig. 4. — Evolution of H ll region with ionizing source off-center from peak density of a 
collapsing region formed from initially uniform gas (model D). Times are given in megayears 
after ionization begins at tjo„. Shown are two-dimensional cuts through the density field 
in the xy-plane containing the peak density. The source is located behind the cut plane. 
Greyscale shows log of density, with values given by the colorbar. Artifacts appear directly 
along grid-lines from the source produced by a slight error in propagating rays in those 
directions. Note the lack of confinement or cometary morphology in the densest regions. 
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Fig. 5. — Time evolution of regions with sources at different positions relative to peak density 
of collapsing region formed from initially uniform gas. Sources are shifted 8 zones (0.125 pc) 
along the {bottom) x-axis, (middle) y and 2;-axes, and (top) x, y, and 2;-axes, (models B, C, 
and D, respectively). Cuts along the xy plane in the plane of the source are shown. (The 
top two panels can be compared to the cuts through the same model in the plane of peak 
density shown in Figure H]). Times are given in megayears, and the greyscale shows log of 
density. Artifacts appear directly along grid-lines from the source produced by a slight error 
in propagating rays in those directions. Morphologies do not depend strongly on position of 
source. 



-33- 




Fig. 6. — Change in peak density over time in models with sources at different positions 
relative to peak density of collapsing region formed from uniform gas. Sources are shifted 8 
zones (0.125 pc) along the x-axis {solid line), y and 2;-axes {dotted line), and x, y, and z-axes 
{dashed line). These are models B, C, and D. 
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Fig. 7. — Column density for the sum of neutral and ionized gas (dominated by the neutral 
gas) in a model of the expansion of an H ll region collapsing from initially uniform gas (model 
D). Projections are shown along the x, y, and 2;-axes as indicated, at the same time as the 
final panel of Figure |H The bright spot is the center of the dense collapsing core, while the 
cross-shaped artifact comes from a slight error in propagating rays directly along grid-lines 
from the source. 
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Fig. 8. — Resolution study of maximum density over time for turbulent models with 64'^ 
{dotted line, model J), 128^ {dashed line, model E), and 256'^ {solid line, model H) zones. 
Equilibrium driven turbulence is first allowed to develop, then gravity is turned on a t a time 
shown by the vertical line on the x-axis. The Jeans criterion (ITruelove et al.l 119971 ) for the 
three different resolutions is shown on the outer y-axis. 
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Fig. 9. — Evolution of H ll region with ionizing source off-center from peak density of a 
collapsing region formed from initially turbulent gas (model H). Times are given in years 
after ionization begins at tjo„. Note change of time units and cube size from Figure |H 
Two-dimensional cuts in the xz plane (for comparison with the xy cuts shown in Figure [TOl) 
centered on the source, with greyscale showing log of density, are shown. 
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Fig. 10. — Resolution study of time evolution of region shown in Figure [H with models E 
and H shown using 128^ and 256^ zones, respectively. Cuts through the xy plane in the plane 
of the source are shown, for comparison in the 256^ case to the xz cut shown in Figure M 
The different resolution models use statistically identical turbulence but not the same actual 
driving pattern, so the shapes should not be compared point by point. Times are again given 
in years, and the greyscale shows log of density. 
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Fig. 11. — Time evolution of regions with sources at different positions relative to peak 
density of collapsing region formed from initially turbulent gas. Sources are shifted 8 zones 
(0.05 pc) along the {bottom) x-axis, {middle) y and 2;-axes, and {top) x, y, and z-axes (models 
E, F, and G, respectively). All models use the same turbulent driving pattern, differing only 
in source position. Cuts through the xy plane in the plane of the source are shown (Note 
that because of the z shift, this is a different plane in the top panels). Times are given in 
years, and the greyscale shows log of density. 
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Fig. 12. — Column density of total gas for a model of the expansion of an H ll region into 
gas collapsing from initially turbulent state (model H). Projections are shown along the x, 
y, and z-axes as indicated, at the same time as the final panel of Figure |9l 
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Fig. 13. — The dimensionless ionizing luminosity C{ri) required to reach a fractional ioniza- 
tion standoff radius ^ = {d — rs)/d for a star a radial distance r — d away from a core with 
density power- law n oc r~'^. The core is to the left, the star to the right, so higher luminosity 
ionizes closer to the center of the core. Solutions are shown for 77 = 3/2 [dashed) and 77 = 2 
[solid). 



